Galton Board
(* Initialization *)
ball = <|
r -> 0.1, (* Radius of the ball *)
n -> 100 (* Number of simultaneous balls *)
|>;
peg = <|
r -> 0.25 (* Radius of the pegs *)
|>;
sim = <|
\[CapitalDelta]t -> 0.05, (* Higher Value need to correct collision positions *)
\[Epsilon] -> 0.4, (* 1.0 means all energy back in the ball *)
g -> {0, -9.81}, (* You can tilt the board *)
n -> 1000 (* Total Number of Balls *)
|>;
(* Initialize Position and Momentum Matrices *)
randomBall := {RandomReal[{-0.05, 0.05}], RandomReal[{0, 1}] + 10};
ball[pos] = Table[randomBall, ball[n]];
ball[mom] = Table[0, ball[n], 2];
ball[bins] = {};
(* Peg position *)
peg[pos] = Flatten[Table[{ii + 0.5 jj, jj}, {ii, -15, 10}, {jj, -1, 8}], 1];
(* Display function *)
positions = ball[pos];
Graphics[{
{LightGray, EdgeForm[Gray], Disk[{#1, #2}, peg[r]]} & @@@ peg[pos],
{PointSize[0.02], Opacity[0.8], Blue, Point[positions // Offload]},
AnimationFrameListener[positions // Offload, "Event"->"frameTick"]
},
AspectRatio -> 0.5,
PlotRange -> {10 {-1, 1}, 10 {0, 1}},
PlotRangeClipping -> True,
ImageSize -> Large
]
EventHandler["frameTick", Function[Null,
If[Length[ball[pos]] > 0,
(* Step *)
(* Scheme: Backward Euler *)
(* Next Step Momentum *)
ball[mom] += Threaded[sim[\[CapitalDelta]t] sim[g]];
(* Next Position Prediction *)
next = ball[pos] + sim[\[CapitalDelta]t] ball[mom];
(* Distance: Do not parallelize, it's already in the algo. *)
dist = DistanceMatrix[next, peg[pos]];
(* Collisions between current and future step *)
col = Parallelize @ MapIndexed[
{v, i} |-> {
SelectFirst[v, (# < ball[r] + peg[r]) & -> "Index"], i[[1]]
} /. ({Missing[_], _} -> Nothing),
dist
];
(* Correct Momentum to include collisions *)
If[Length[col] > 1,
{
up, val
} = Transpose[
Parallelize @ MapApply[
{p, b} |-> {
b,
sim[\[Epsilon]] * Norm[ball[mom][[b]]] Normalize[next[[b]] - peg[pos][[p]]]
},
col
]
];
momt = ball[mom];
momt[[up]] = val;
ball[mom] = momt;
];
(* Update pos with corrected momentum *)
ball[pos] += sim[\[CapitalDelta]t] ball[mom];
(* Filter through the balls, remove/reset and register those exiting *)
exited = Select[ball[pos], #[[2]] < 0 & -> "Index"];
ingame = Select[ball[pos], #[[2]] > 0 & -> "Index"];
If[Length[exited] > 0,
ball[bins] = Join[ball[bins], ball[pos][[exited, 1]]]
];
(* Remove them from the pos and mom matrices *)
ball[pos] = ball[pos][[ingame]];
ball[mom] = ball[mom][[ingame]];
(* Insert new ball *)
left = sim[n] - Length[ball[bins]] - Length[ball[pos]];
toAdd = Min[left, ball[n] - Length[ball[pos]]];
Do[
AppendTo[ball[pos], randomBall];
AppendTo[ball[mom], {0, 0}],
toAdd
];
positions = ball[pos];
]]];